Cluster Analysis
# Hierarchical Clustering based on square root transformed abundance data
# Agglomeration used group average (= UPGMA) method
d <- vegdist(b^0.25)
hc <- hclust(d, method="average")
#plot(x = hc, labels = row.names(hc), cex = 0.5, hang=-1)
#rect.hclust(tree = hc, k = 4)
fac <- ta[, 1:4]
fac$Group <- cutree(hc, 4)
# Reorder group from north to south as Groups 1, 2, 3
fac$Group <- factor(fac$Group, levels=c(1,2,3,4), labels = c(2, 3, 1, 4))
# Convert dendrogram to ggplot style
dhc <- as.dendrogram(hc)
ghc <- dendro_data(dhc, type="rectangle")
# Merge sample info to dendrogram label data frame
ord <- match(label(ghc)$label, rownames(b))
ghc[["labels"]] <- cbind(ghc[["labels"]], fac[ord,])
# Group membership
# Group 1 (Shitiping = 82.6%)
table(subset(fac, Group==1)$Site)
##
## STP JH JLL
## 19 2 2
# Group 2 (Jihuei = 52.6%, Jialulan = 28.9%)
table(subset(fac, Group==2)$Site)
##
## STP JH JLL
## 7 20 11
# Group 3 (Jialulan = 77.6%, Jihuei = 17.2%)
table(subset(fac, Group==3)$Site)
##
## STP JH JLL
## 3 10 45
Nonmetric Multidimensional Scaling
# Remove outlier JH_AP_Jania sp.2_1
md <- metaMDS(vegdist(b[-23,]^0.25))
## Run 0 stress 0.1395381
## Run 1 stress 0.1493668
## Run 2 stress 0.1431431
## Run 3 stress 0.1454937
## Run 4 stress 0.142695
## Run 5 stress 0.1418538
## Run 6 stress 0.1418331
## Run 7 stress 0.1380858
## ... New best solution
## ... Procrustes: rmse 0.04631679 max resid 0.1232183
## Run 8 stress 0.1383089
## ... Procrustes: rmse 0.01505367 max resid 0.07531758
## Run 9 stress 0.1397265
## Run 10 stress 0.14193
## Run 11 stress 0.1406395
## Run 12 stress 0.1397602
## Run 13 stress 0.1382344
## ... Procrustes: rmse 0.004775539 max resid 0.04929032
## Run 14 stress 0.1383071
## ... Procrustes: rmse 0.01515002 max resid 0.07518309
## Run 15 stress 0.1398386
## Run 16 stress 0.1380858
## ... New best solution
## ... Procrustes: rmse 3.250285e-05 max resid 0.0002187851
## ... Similar to previous best
## Run 17 stress 0.1395381
## Run 18 stress 0.1414666
## Run 19 stress 0.1389654
## Run 20 stress 0.1425775
## *** Best solution repeated 1 times
stress <- paste("Stress = ", deparse(round(md$stress,2)))
keep <- names(colSums(b[-23,])%>%sort(decreasing = TRUE))[1:5]
blab <- wascores(md$points, b[-23,]^0.25, expand = TRUE)[keep,] %>% as.data.frame
blab$label <- paste("italic(", sub(" ", "~~", rownames(blab)), ")", sep="")
# Only fits environmental vector to the replicate with grain size data
keep <- !env[-23, 1:4]%>%rowSums%>%is.na
set.seed(100)
(fit <- envfit(md$points[keep,], env[-23,][keep,]))
##
## ***VECTORS
##
## MDS1 MDS2 r2 Pr(>r)
## Silt2Clay 0.29515 0.95545 0.0264 0.475
## Fine2VeryFine 0.53410 0.84542 0.1250 0.031 *
## Medium 1.00000 0.00267 0.2581 0.001 ***
## Coarse2VeryCoarse 0.88051 -0.47403 0.0912 0.078 .
## Temperature -0.69724 0.71683 0.0690 0.158
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
evec <- as.data.frame(scores(fit, display = "vectors")) * ordiArrowMul(fit, fill = 1.5)
evec$label <- row.names(evec)
# Mantel tests on each variables
man <- list()
man[[1]] <- mantel(vegdist(b[-23,][keep,]^0.25), dist(env[-23,][keep, "Silt2Clay"]))
man[[2]] <- mantel(vegdist(b[-23,][keep,]^0.25), dist(env[-23,][keep, "Fine2VeryFine"]))
man[[3]] <- mantel(vegdist(b[-23,][keep,]^0.25), dist(env[-23,][keep, "Medium"]))
man[[4]] <- mantel(vegdist(b[-23,][keep,]^0.25), dist(env[-23,][keep, "Coarse2VeryCoarse"]))
man[[5]] <- mantel(vegdist(b[-23,]^0.25), dist(env[-23, "Temperature"]))
evec$mantel_r <- lapply(man, FUN=function(x)x$statistic) %>% unlist
evec$signif <- lapply(man, FUN=function(x)x$signif) %>% unlist
#evec$pvals <- cut(evec$signif, breaks=c(0, 0.01, 0.05, 0.1, 1), labels=c("<0.01", "<0.05", "<0.1", "<1")) %>% factor
evec$pvals <- factor(evec$signif)
print(evec)
## MDS1 MDS2 label mantel_r signif
## Silt2Clay 0.1415600 0.458253191 Silt2Clay -0.05470974 0.827
## Fine2VeryFine 0.5575815 0.882594237 Fine2VeryFine 0.07029091 0.083
## Medium 1.5000000 0.003998335 Medium 0.21036458 0.001
## Coarse2VeryCoarse 0.7850268 -0.422627016 Coarse2VeryCoarse 0.11119698 0.006
## Temperature -0.5406990 0.555889465 Temperature 0.13973555 0.001
## pvals
## Silt2Clay 0.827
## Fine2VeryFine 0.083
## Medium 0.001
## Coarse2VeryCoarse 0.006
## Temperature 0.001
# Merge MDS to environmental data frame
p2 <- ggplot(data=cbind(md$points, fac[-23,]), aes(x=MDS1, y=MDS2))+
geom_point(alpha=1, size=5, stroke=1, aes(colour=Site, shape=Season))+
geom_mark_hull(concavity = 10, aes(group=Group), colour="black", linetype=2)+
geom_segment(data= evec, aes(x=0, y=0, xend=MDS1, yend=MDS2, linewidth=pvals, linetype=pvals),
arrow = arrow(length=unit(.4, 'cm')), colour="red")+
annotate("text", x=-0.9, y=1.1, label=stress, size=5) +
annotate("text", x=1, y=0, label="1", size=10) +
annotate("text", x=-0.5, y=-0.6, label="2", size=10) +
annotate("text", x=-0.2, y=0.5, label="3", size=10) +
annotate("text", x=1.7, y=1.2, label="b", size=7) +
scale_shape_manual(values=c(1, 19))+
scale_color_viridis_d()+
scale_linewidth_manual(values=c(2, 1, 0.5, 0.5))+
scale_linetype_manual(values=c(1, 1, 1, 2))+
geom_label_repel(data=blab, aes(x=MDS1, y=MDS2, label=label), colour="black",
fill=gray(1, 0.6), size=4, label.padding = unit(0.25, "lines"), parse=T)+
geom_label_repel(data=evec, aes(x=MDS1*1.2, y=MDS2*1.2, label=label), colour="red",
fill=gray(1, 0.6), size=4, label.padding = unit(0.25, "lines"), force=2)+
labs(linewidth="p-value", linetype="p-value")+
theme_bw() %+replace% large
ggdendrogram(hc, rotate = FALSE)+
geom_point(data=label(ghc), aes(x, y, colour=Site, shape=Season), size=2, stroke=0.5)+
geom_hline(yintercept = 0.8, linetype=2, colour="red")+
annotate("text", x = 95, y = 0.85, label = "2", size = 5)+
annotate("text", x = 42, y = 0.85, label = "3", size = 5)+
annotate("text", x = 10, y = 0.85, label = "1", size = 5)+
scale_shape_manual(values=c(1, 19))+
scale_color_viridis_d()+
labs(x= "Replicate",y = "Bray-Curtis Dissimilarity")+
theme_bw() %+replace% theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), axis.line.y = element_line(colour="black")) %+replace% rotate
m <- acast(dat=ta, Microhabitat~Site+Season, fun.aggregate = length)
hc <- vegdist (m) %>% hclust(method="average")
clust_names <- row.names(m)[rev(hc$order)]
out <- cbind(md$points, fac[-23,])
out$Microhabitat <- factor(out$Microhabitat, levels=clust_names)
# Color code the microhabitat
cols <- c("#222222", "#f3c300", "#875692", "#f38400", "#a1caf1", "#be0032", rep("#c2b280", 7), rep("#848482", 3),"#008856")
shapes <- c(rep(19, 6), 0:6, 15, 17, 18:19)
# MDS for each site and season
md <- lapply(splitBy(~Site+Season, ta), FUN=function(x)metaMDS(x[,-1:-4]))
## Wisconsin double standardization
## Run 0 stress 0.09177207
## Run 1 stress 0.09177187
## ... New best solution
## ... Procrustes: rmse 0.000338808 max resid 0.001136196
## ... Similar to previous best
## Run 2 stress 0.09240116
## Run 3 stress 0.1102589
## Run 4 stress 0.09185124
## ... Procrustes: rmse 0.07012022 max resid 0.1591274
## Run 5 stress 0.09901149
## Run 6 stress 0.09185037
## ... Procrustes: rmse 0.06994872 max resid 0.1586776
## Run 7 stress 0.1204544
## Run 8 stress 0.09857996
## Run 9 stress 0.09058073
## ... New best solution
## ... Procrustes: rmse 0.1173281 max resid 0.306095
## Run 10 stress 0.09240098
## Run 11 stress 0.1198589
## Run 12 stress 0.09185052
## Run 13 stress 0.0955431
## Run 14 stress 0.09050981
## ... New best solution
## ... Procrustes: rmse 0.08364446 max resid 0.2192915
## Run 15 stress 0.0909706
## ... Procrustes: rmse 0.0747375 max resid 0.2755401
## Run 16 stress 0.1253892
## Run 17 stress 0.125494
## Run 18 stress 0.09050979
## ... New best solution
## ... Procrustes: rmse 0.0006690435 max resid 0.001675329
## ... Similar to previous best
## Run 19 stress 0.09058075
## ... Procrustes: rmse 0.08327585 max resid 0.2200089
## Run 20 stress 0.09058272
## ... Procrustes: rmse 0.08290963 max resid 0.2208044
## *** Best solution repeated 1 times
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1345332
## Run 1 stress 0.1480562
## Run 2 stress 0.1468133
## Run 3 stress 0.1551481
## Run 4 stress 0.1383558
## Run 5 stress 0.153372
## Run 6 stress 0.1322342
## ... New best solution
## ... Procrustes: rmse 0.09070836 max resid 0.298377
## Run 7 stress 0.1403992
## Run 8 stress 0.1755181
## Run 9 stress 0.1678867
## Run 10 stress 0.1564364
## Run 11 stress 0.1442084
## Run 12 stress 0.1524216
## Run 13 stress 0.1361397
## Run 14 stress 0.1403993
## Run 15 stress 0.17642
## Run 16 stress 0.153372
## Run 17 stress 0.1346216
## Run 18 stress 0.1678868
## Run 19 stress 0.1550996
## Run 20 stress 0.139706
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 20: stress ratio > sratmax
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1208441
## Run 1 stress 0.09427342
## ... New best solution
## ... Procrustes: rmse 0.1234327 max resid 0.2720241
## Run 2 stress 0.1098328
## Run 3 stress 0.09427323
## ... New best solution
## ... Procrustes: rmse 7.529209e-05 max resid 0.0002953634
## ... Similar to previous best
## Run 4 stress 0.1382285
## Run 5 stress 0.1349709
## Run 6 stress 0.09332658
## ... New best solution
## ... Procrustes: rmse 0.03245163 max resid 0.1489107
## Run 7 stress 0.1196955
## Run 8 stress 0.09427334
## Run 9 stress 0.1212999
## Run 10 stress 0.125197
## Run 11 stress 0.1175232
## Run 12 stress 0.1265351
## Run 13 stress 0.1170811
## Run 14 stress 0.09332658
## ... New best solution
## ... Procrustes: rmse 0.000273151 max resid 0.0009519536
## ... Similar to previous best
## Run 15 stress 0.1099111
## Run 16 stress 0.09427323
## Run 17 stress 0.09332666
## ... Procrustes: rmse 0.0003160555 max resid 0.001152594
## ... Similar to previous best
## Run 18 stress 0.1145061
## Run 19 stress 0.1054994
## Run 20 stress 0.1294716
## *** Best solution repeated 2 times
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.138881
## Run 1 stress 0.1380855
## ... New best solution
## ... Procrustes: rmse 0.02290013 max resid 0.09354521
## Run 2 stress 0.2044978
## Run 3 stress 0.2363496
## Run 4 stress 0.1750268
## Run 5 stress 0.2121027
## Run 6 stress 0.2069126
## Run 7 stress 0.1518363
## Run 8 stress 0.2125752
## Run 9 stress 0.1393694
## Run 10 stress 0.216466
## Run 11 stress 0.1553091
## Run 12 stress 0.133009
## ... New best solution
## ... Procrustes: rmse 0.08724517 max resid 0.3626722
## Run 13 stress 0.1638621
## Run 14 stress 0.1400827
## Run 15 stress 0.2068781
## Run 16 stress 0.2026245
## Run 17 stress 0.1564282
## Run 18 stress 0.1526192
## Run 19 stress 0.2163723
## Run 20 stress 0.2174266
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 15: stress ratio > sratmax
## 5: scale factor of the gradient < sfgrmin
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.05720906
## Run 1 stress 0.0572091
## ... Procrustes: rmse 0.0002015069 max resid 0.0003602265
## ... Similar to previous best
## Run 2 stress 0.0572092
## ... Procrustes: rmse 0.0002634092 max resid 0.0004780251
## ... Similar to previous best
## Run 3 stress 0.05720917
## ... Procrustes: rmse 0.0002193702 max resid 0.0004674462
## ... Similar to previous best
## Run 4 stress 0.05720909
## ... Procrustes: rmse 7.852231e-05 max resid 0.00014762
## ... Similar to previous best
## Run 5 stress 0.05720917
## ... Procrustes: rmse 0.0002701912 max resid 0.0004919481
## ... Similar to previous best
## Run 6 stress 0.05720914
## ... Procrustes: rmse 0.0002475757 max resid 0.0004502463
## ... Similar to previous best
## Run 7 stress 0.05720912
## ... Procrustes: rmse 0.0002307228 max resid 0.0004188804
## ... Similar to previous best
## Run 8 stress 0.08766255
## Run 9 stress 0.08093205
## Run 10 stress 0.08093206
## Run 11 stress 0.1086705
## Run 12 stress 0.05720904
## ... New best solution
## ... Procrustes: rmse 4.555359e-05 max resid 7.944265e-05
## ... Similar to previous best
## Run 13 stress 0.05720907
## ... Procrustes: rmse 0.0001330564 max resid 0.0002378443
## ... Similar to previous best
## Run 14 stress 0.08093209
## Run 15 stress 0.05720903
## ... New best solution
## ... Procrustes: rmse 3.619008e-05 max resid 0.0001015064
## ... Similar to previous best
## Run 16 stress 0.05720904
## ... Procrustes: rmse 2.007738e-05 max resid 5.202092e-05
## ... Similar to previous best
## Run 17 stress 0.07810263
## Run 18 stress 0.08093205
## Run 19 stress 0.08093206
## Run 20 stress 0.05720908
## ... Procrustes: rmse 0.0001108149 max resid 0.0002103537
## ... Similar to previous best
## *** Best solution repeated 3 times
## Wisconsin double standardization
## Run 0 stress 6.522281e-33
## Run 1 stress 0
## ... New best solution
## ... Procrustes: rmse 0.1307195 max resid 0.2126991
## Run 2 stress 0
## ... Procrustes: rmse 0.1368719 max resid 0.2230256
## Run 3 stress 0
## ... Procrustes: rmse 0.1450721 max resid 0.287975
## Run 4 stress 0
## ... Procrustes: rmse 0.1626683 max resid 0.2638379
## Run 5 stress 0
## ... Procrustes: rmse 0.09148311 max resid 0.1643084
## Run 6 stress 1.910875e-05
## ... Procrustes: rmse 0.149761 max resid 0.2413231
## Run 7 stress 0
## ... Procrustes: rmse 0.1231937 max resid 0.2250026
## Run 8 stress 0
## ... Procrustes: rmse 0.1865359 max resid 0.3227937
## Run 9 stress 0
## ... Procrustes: rmse 0.1331486 max resid 0.2704025
## Run 10 stress 0
## ... Procrustes: rmse 0.115681 max resid 0.1553413
## Run 11 stress 0
## ... Procrustes: rmse 0.1302713 max resid 0.2177031
## Run 12 stress 6.264849e-05
## ... Procrustes: rmse 0.1032593 max resid 0.1609292
## Run 13 stress 0
## ... Procrustes: rmse 0.139783 max resid 0.2708793
## Run 14 stress 0
## ... Procrustes: rmse 0.1839108 max resid 0.3134858
## Run 15 stress 0
## ... Procrustes: rmse 0.1270628 max resid 0.1964893
## Run 16 stress 0
## ... Procrustes: rmse 0.172901 max resid 0.2816818
## Run 17 stress 0
## ... Procrustes: rmse 0.1588443 max resid 0.2648712
## Run 18 stress 0
## ... Procrustes: rmse 0.1501897 max resid 0.2424697
## Run 19 stress 0
## ... Procrustes: rmse 0.1613569 max resid 0.2847132
## Run 20 stress 0
## ... Procrustes: rmse 0.1396997 max resid 0.2310131
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 20: stress < smin
# Stress levels
st <- lapply(md, FUN=function(x)x$stress) %>% ldply
st <- cbind(strsplit(st$.id, split="[|]") %>% ldply, st[, -1])
names(st) <- c("Site", "Season", "stress")
st$stress <- paste("Stress = ", round(st$stress,2))
st$Site <- factor(st$Site, levels=c("STP", "JH", "JLL"))
# MDS scores
md <- cbind(ta[,1:4], lapply(md, FUN=function(x)x$points) %>% ldply)
md$Microhabitat <- factor(md$Microhabitat, levels = clust_names)
ggplot(data=md, aes(x=MDS1, y=MDS2, colour=Microhabitat, shape=Microhabitat))+
geom_point(size=3)+
geom_text(data = st, inherit.aes = FALSE, aes(x=Inf, y = Inf, label = stress), vjust=1.2, hjust=1.05)+
scale_colour_manual(values=cols)+
scale_shape_manual(values=shapes)+
facet_grid(Site~Season, scales="free")+
theme_bw() %+replace% large
